%be sure you use the saved variables in the sequence!!!
%inputs: location, units of t_tof , units t_hold
%units are strings: 's', 'ms', 'us'
%defaults are ms
    % Fits a double Gaussian function to data for finding the temperature of
    % atoms released from the dipole trap. Also attempts to measure the
    % temperature of the released MOT cloud for comparison.
    % if after tHOLD the MOT falls a distance larger that the 
    %imaged space, only one Gaussian is fitted






function  fit_tof = MeasureTemperatureTOF(location,t_tof_scale,t_hold_scale)

    

    %location='13_02_55';

    expp = ParamImgSys(GetImgConfig(location));

        switch t_tof_scale
            case 's'
                t_tof_scale=1;
            case 'ms'
                t_tof_scale=1e-3;
            case 'us'
                t_tof_scale=1e-6;
        otherwise
            error('Check units');
                    
        end
        
         switch t_hold_scale
            case 's'
                t_hold_scale=1;
            case 'ms'
                t_hold_scale=1e-3;
            case 'us'
                t_hold_scale=1e-6;
        otherwise
            error('Check units');
                    
        end
    %t_tof_scale = 1e-6; %scale for t_tof. Saved value is not in mks units necessarily 
    
    
    %t_hold_scale = 1e-3; %scale factor for t_hold
    position_scale = expp.px_size/expp.mag;


    %ROIx =  538:1013; 
    %

    numRuns = length(dir([location,'/*.h5']));
    cons = Constants();
    yC = 1:expp.short_px_nb;
    xC = 1:expp.long_px_nb;

    getData = 1;
    t = [];
    ODav = [];
    if getData
        for n = 2:numRuns
            n
            tN = ReadVariable(location , n-1 , 't_tof' ) * t_tof_scale; %in s

            %tN = (double(h5read([path '\run_' num2str(n-1) '.h5'],'/VariableValues/t_tof')))*1e-6;
            tHOLD = ReadVariable(location , n-1 , 't_hold' ) * t_hold_scale; %in s
            %tHOLD = (double(h5read([path '\run_' num2str(n-1) '.h5'],'/VariableValues/t_hold')))*1e-3;
            new_t = find(t == tN);
            ODint = ReadDataAnalysisRes( location , n-1,'xData');
            %ODint = double(h5read([path '\run_' num2str(n-1) '.h5'],'/dataAnalysis/xData'));
            if length(new_t) == 0 % if it's a 'new' time, cat. to data
                new_t = length(t) + 1;
                t = [t tN];
                ODav = [ODav; ODint];
            else % else average with existing data
                ODav(new_t,:) = ODav(new_t,:) + ODint;
            end
        end
        save vars
    else
        load vars
    end

    %initial guesses

    x = [0:length(ODav)-1]*position_scale;
    x_mot_m_ig = tHOLD^2*cons.g0/2; %in m
    x_mot_pix_ig = round(tHOLD^2*cons.g0/position_scale); %in pixels

    %assuming thes values for initial guess on the MOT signal
    sigma_mot_ig=0.5e-3;
    T_mot_ig=40e-6;%K MOT temp
    %same for atoms in ODT
    sigma_dip_ig=50e-6; %m



    for I = 1:length(t)
        sig = ODav(I,:)';

        %fit two Gaussians 
        %initial guess: initial values for dipole gaussian are given in 
        %terms of the MOT temperature and falling time (tHOLD)
        %if after tHOLD the MOT falls a distance larger that the 
        %image view, only one Gaussian is fitted


        [aux,aux_index] = max(sig);

        p_ig(1) = aux;
        p_ig(2) = x(aux_index);
        p_ig(3) = sigma_dip_ig;

        aux2 = round(p_ig(2)-x_mot_pix_ig);

        options = optimset('Display','off','TolFun',1e-16,'TolX',1e-16);

        if aux2 > 0 %fit double Gaussian

            p_ig(4) = sig(aux2);
            p_ig(5) = p_ig(2)-x_mot_pix_ig;
            p_ig(6) = sigma_mot_ig;
            p_ig(7) = min(sig);
            fit_x(I,:) = fminsearch(@(p)FitDoubleGauss(p,x,sig), ...
             [p_ig(1:7)],options);

        else %fit single Gaussian
            p_ig(4) = (sig(1)-sig(end))/2;
            fit_x(I,:) = fminsearch(@(p)FitGauss(p,x,sig), ...
                [p_ig(1:4)],options);

        end




    end

    figure(1)
    %plot(x,ODav')
    hold on

    sigma_x_dip=[];
    peak_gaussian_dip=[];
    sigma_x_mot=[];
    peak_gaussian_dip=[];
    for I = 1:length(t)
        params = fit_x(I,:);
        if length(params)>4
            consideringmot = 1;
            plot(x,DoubleGauss(params,x),x,ODav(I,:)')
            sigma_x_mot=[sigma_x_mot, param(6)];
            peak_gaussian_mot=[peak_gaussian_mot,params(4)];
        else
            consideringmot = 0;
            plot(x,gauss(params,x),x,ODav(I,:)');
        end
        sigma_x_dip=[sigma_x_dip,params(3)]
        peak_gaussian_dip=[peak_gaussian_dip,params(2)];
    end


    hold off
    legend(num2str(t'))
    xlabel('x (m)')



   % temperature of atoms in ODT
    for count = 1:length(sigma_x_dip)
        if sigma_x_dip(count) > 2.6e-4
            sigma_x_dip(count) = 0;
            peak_gaussian_dip(count) =0;
            t(count)=0;
        end
    end
     sigma_x_dip= sigma_x_dip( sigma_x_dip~=0);
     peak_gaussian_dip=peak_gaussian_dip(peak_gaussian_dip~=0);
     t=t(t~=0);
     
    fit_tof_coeff_dip = polyfit(t.^2,sigma_x_dip.^2,1);
    sigma2_v_dip = fit_tof_coeff_dip(1);
    sigma2_x0_dip = fit_tof_coeff_dip(2);
    figure(2)
    plot(t.^2, sigma_x_dip.^2, 'b.' , t.^2,t.^2*sigma2_v_dip+sigma2_x0_dip , 'g-')
    title('Temperature measurement')
    xlabel('t^2 (s^2)')
    ylabel('\sigma _x ^2 (m^2)')
    
    T_fit_dip = cons.massRb87kg/cons.kB*sigma2_v_dip

    %calculating g
    fit_g_coeff_dip =  polyfit(t.^2,peak_gaussian_dip,1);
    g_fit_dip = fit_g_coeff_dip(1)*2

    %plot(t.^2,peak_gaussian_dip,t.^2,t.^2*fit_g_coeff_dip(1)+fit_g_coeff_dip(2))

    if consideringmot
        t=t+tHOLD;
        fit_tof_coeff_mot = polyfit(t.^2,sigma_x_mot.^2,1);
        sigma2_v_mot = fit_tof_coeff_mot(1);
        sigma2_x0_mot = fit_tof_coeff_mot(2);

        %plot(t.^2, sigma_x_mot.^2,t.^2,t.^2*sigma2_v_mot+sigma2_x0_mot)

        T_fit_mot = cons.massRb87kg/cons.kB*sigma2_v_mot

        %calculating g
        fit_g_coeff_mot =  polyfit(t.^2,peak_gaussian_mot,1);
        g_fit_mot = fit_g_coeff_mot(1)*2

        %plot(t.^2,peak_gaussian_mot,t.^2,t.^2*fit_g_coeff_mot(1)+fit_g_coeff_mot(2))


    end